;;---------------------------------------------------------------------  
;; creat zonal mean radon and ipr distribution 
;;---------------------------------------------------------------------  
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" 
load "$KAILIBROOT/lib_interp.ncl"
load "$KAILIBROOT/lib_radon.ncl"


begin
  
  overland = True 

  path = "/pf/m/m222044/libdata/"
  pata = "/pf/m/m222044/lib/ionprod/figure_geograph/input/"

  season = "JJA" 
  run = "R46H" 

  load "load_interp" 

  filenama = path+"mbase_T63L31.nc" 
  filenamb = pata+run+"_IPRC_"+season+".nc" 
  filenamc = path+"geop_T63L31_"+season+".nc"   

  print("reading "+filenama) 
  print("reading "+filenamb) 
  print("reading "+filenamc) 

  fila = addfile(filenama,"r") 
  filb = addfile(filenamb,"r") 
  filc = addfile(filenamc,"r")  

  time     = fila->time
  lat      = fila->lat
  lon      = fila->lon
  phis     = fila->geosp
  slm      = fila->slm
  hyam     = filb->hyam
  hybm     = filb->hybm
  rns      = filb->IPR
  geop     = filc->var156 
 
  print("--- data read ---") 
  


  ;; new height levels 

  hh = (/4000.,3500.,3000.,2750.,2500.,2250.,2000.,1750.,1500.,1250.,1000.,900.,800.,700.,600.,500.,400.,300.,250.,200.,150.,100.,60.,30.,0./) 
  
  nlat = dimsizes(lat) 
  nlon = dimsizes(lon) 
  ntim = dimsizes(time) 
  nlev = dimsizes(hh) 

  rnh = new((/ntim,nlev,nlat,nlon/),"float") 
  rnh = 0. 
  rnh@_FillValue = -999. 

  rnh!0 = "time" 
  rnh!1 = "lev" 
  rnh!2 = "lat" 
  rnh!3 = "lon" 

  rnh&time= time
  rnh&lev = hh
  rnh&lat = lat
  rnh&lon = lon

  rmh = rnh
  roh = rnh
  iph = rnh

  print("--- array created ---") 

  phisp = conform(geop,geop(:,30,:,:),(/0,2,3/)) 
 ;phisp = conform(geop,phis(:,:,:),   (/0,2,3/)) 

  ;; height above ground (m)  

  geop = geop - phisp + 60. ;;/ 10. 

  print("max,min of geop : "+max(geop)+"   "+min(geop)) 
  print("--- height caled ---") 

  ;;do m = 0, ntim - 1 

  ;;print("m = "+m) 

  ;; ion production rate 
 


  do j = 0, nlat-1
  do i = 0, nlon-1
      rnh(:,:,j,i) = (/ linint1( geop(:,:,j,i), rns(:,:,j,i), False, hh, 0 ) /)  
  end do
  end do

  do k = 0, nlev-1  
     rnh(:,k,:,:) = where(ismissing(rnh(:,k,:,:)),rns(:,30,:,:),rnh(:,k,:,:))    
  end do 

 
  ;;end do


  if (overland) then 

  ;do j = 0, nlat - 1
  ;   if(all(slm(:,j,:).le.1.e-5))  then 
  ;      slm(:,j,:) = 1. 
  ;   end if  
  ;end do

  slmp = conform(rnh,slm,(/0,2,3/)) 

  ;printVarSummary(slmp) 
  ;printVarSummary(rnh) 

  rnh = where(slmp.lt.0.5,-999.,rnh)  

  end if ;; overland 


  print("--- output ---") 


;; #/m3/s -> #/cm3/s 

  rnh = rnh * 1.E-06 
   

 ;filenamo = "out_all_"+run+"_"+season+".nc" 
  filenamo = "out_IPRC_"+run+"_"+season+".nc" 

  system("rm "+filenamo) 
  
  filo = addfile(filenamo,"c") 

  filo->IPRC=rnh  
  filo->geop=geop 

end



